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Abstract Modelling isolated rotating stars at any rotation rate is a challenge 
for the next generation of stellar models. These models will couple dynamical 
aspects of rotating stars, like angular momentum and chemicals transport, 
with classical chemical evolution, gravitational contraction or mass-loss. Such 
modelling needs to be achieved in two dimensions, combining the calculation 
of the structure of the star, its mean flows and the time-evolution of the whole. 
We present here a first step in this challenging programme. It leads to the 
first self-consistent two-dimensional models of rotating stars in a steady state 
generated by the ESTER code. In these models the structure (pressure, den- 
sity and temperature) and the flow fields are computed in a self-consistent 
way allowing the prediction of the differential rotation and the associated 
meridian circulation of the stars. After a presentation of the physical proper- 
ties of such models and the numerical methods at work, we give the first grid 
of such models describing massive and intermediate-mass stars for a selection 
of rotation rates up to 90% of the breakup angular velocity. 
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1 Introduction 



1.1 The astrophysical context 

In the last ten years rotation has become an essential part of stellar models. 
It influences all stages of stellar evolution. For instance, during the formation 
process, a drastic reduction of the specific angular momentum, by a factor of 
order 10 5 , is observed when matter passes from the initial molecular cloud to 
the newly-born main-sequence star. Later, during the evolution of the star 
on the main sequence, many hydrodynamical instabilities, along with merid- 
ian circulations, drive some mixing in the radiative zones. The effects of this 
mixing are actually observed on the surface of many stars, which show ele- 
ments obviously synthesi zed in their core region. In addition, recent work by 



Frischknecht et al.l (|2012l ) shows that rotation plays a role in the nucleosyn- 
thesis of s-elements in massive stars. Effects of rotation are also thought to be 
important in understanding the evolution and yields of the first generation 
of stars. Lacking in metals, these stars were more compact than present stars 
and had also weaker winds. Thus, their rotation was certainly faster than 
that of present day stars. It is therefore crucial to master rotational effects 
in order to reconstruct the history of metal enrichment in galaxies and to 
understand how the observed stars have been influenced by the first genera- 
tion. Lastly, we may mention an example of the importance of rotation in the 
end of the life of stars: recent works on the gravitational collapse of massive 



stars (e.g. iMetzger et al.M201 If) insist on the importance of the combined ef- 
fect of rotation and magnetic fields to model the explosion of supernovae and 
the associated production of a gamma-ray burst. The few foregoing examples 
are just illustrative, since rotation influences many other aspects of stellar 
physics like the oscillation spectrum, mass-loss etc. 



1.2 The ID-models 



Presently, rotation is inclu ded in spher ically symmetric stellar models through 
the approach proposed by|Zahn|(l992). Since these models are one-dimensio- 



nal all of the effects of rotation are averaged on spheres. Thus differential 
rotation is only in the radial direction (said to be shellular) and no meridian 
flow appears explicitly. Because only the first terms of the spherical harmonic- 
expansion of fields are included, this modelling is valid at small rotation rates. 
An equally important part of the models is that they assume the existence 
of some small-scale turbulence in the radiative regions that efficiently trans- 
port momentum horizontally (i.e. on isobars), erasing latitudinal gradients of 
angular velocity. 
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Although limited to slow rotation and needing some ad hoc constants (like 
turbulent diffusivities) , this approach has the great merit of allowing the use 
of ID-models, which are now very precise as far as microphysics is concerned. 
It thus allowed investigators to reproduce many obser ved features of stars: 
surfa ce abundance of lithium as a function mass (e.g. ICharbonnel fe; Talon 
1999h . t he (relative) high number of red super-giants in low-metallicity galax- 



ies (e.g. iMaeder fc Mevnetll200ll), o r the ratio of type Ibc to type II super 
novae re.g. lMevnetfcMaederil2005h . etc. 

Although these examples reflect a truly successful modelling, the under- 
standing of the effects of rotation is still incomplete because in many circum- 
stances rotation is fast. Stellar conditions thus do not meet the requirements 
of the models and the use of Zahn's approach becomes dubious. Clearly, we 
are missing a more detailed view of reality. For instance, no one knows the 
rotation rates that are allowed in the foregoing ID models. 



1 . 3 The history of 2D-models of rotating stars 



The solution to the deficiencies of ID-models about rotation will come from 
the use of multi-dimensional models (two-dimensional at least), which in- 
clude the mean-field hydrodynamics along with the centrifugal distortion of 
the star. Unfortunately, building such models turned out to be a very difficult 
challenge. The story of this modelling is illuminatin g and we briefly su mma- 
rize it below, following the more detailed account of Rieutordl (l2006bl). 

The first step dates back to polytropic hydrostatic models of lJamesl (|l964l) . 
An important st ep forward was made by a US group around P. Bodenheimer 
and J. Ostriker (lOstriker fe Marklll968tlOstriker fc Bodenheimerlll968tlMark 
1968HOstriker fc Hartwickll968HJacksonll970l:lBodenheimer fc Ostrikerll970 ; 



Bodenheimerl 1 1 9 7 It iBodenheimer fc Ostrikerl I1973T ) . Their works are based 



on the Self- Consistent Field method which solves Poisson's equation for the 
gravitational potential, A<p = AnGp, with the Green function, namely 



-G 



P(x') 



dh 



This solution readily includes the boundary conditions at infinity for the 
gravitational field. 

Making general stellar models was impeded by many numerical difficulties. 
The codes could not deal with stellar masses less than 9 M nor with very 
rapid rotation. An indication of the precision reached by these models is given 
by the virial test (see below), fjackson ( 1970l ) found it to be around 4xl0 -3 . 

Soon after, M. Clement took up the challenge and solved directly the Pois- 
son equation with a finite difference discretization (|Clementlll974l ). Results 
improved as testified by a better virial test at 2xlCP 4 . 
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Another series of work aimed at the construction of 2D-models was 



published by the Japanes e school led by Y. Eriguchi (see lEriguchi 1978 
Eriguchi fc Sugimotolll981 ). These wor ks led to the first attempts to account 



for the baroclinicity of radiative zones (jUrvu fc Eriguchilll994 119951 ). but the 
neglecting of viscosity imposed the prescription o f some ad hoc constraints. 
The series ended with the work of IShindo et al.l ( 19971 ) who improved the 



microphysics but relaxed the dynamics, considerin g only barotropic models. 

To be complete, we need to mention the work of lEriguchi fc Mullen ([1985, 
19911) who attacked the problem in a different way, using a mapping of the 
star. The grid points are indeed mapped to the sphere through the relation 



n(o k ) = bR.(o k ) 

where £ is a new radial variable spanning [0,1]. Finite differences are used 
for both variables ( and 8. The ensuing code was quite robust, being able to 
compute configurations quite far from the pure sphere, but the precision of 
calculation as monitored by the virial test remained around 5 x 10~ 4 . 

The more recent effort s on 2D modell ing not related to the ESTER 
project have been those of Roxburgh! ( 2004 ). motivated by th e need of mod 



els of rapidl y rotating stars for asteroseismology, and those of Uackson et al 



(|2004l 120051) motiva ted by the discovery of the extre me flattening of the Be 



star Achernar (e.g. iDomiciano de Souza et al.l 120051 ) . These two works use 



barotr opic models with an imposed rotation law. We note that Uackson et al 



(2004) used a new version of the Self-Consistent Field technique which leads 
to a much more robust code, not restricted to a given mass range. These 
models have also been used t o investigate the ac oustic osci llation sp e ctrum 
of rapidly rotating stars by ( Reese et al. 2009bl iah. Lastly, iDeupreej ( 2011 ) 
worked out a series of 2D-models with uniform rotation for masses between 
1.625 and 8 M©. 



2 The route to ideal models 



As may be observed, the weakest point of the foregoing models is their lack of 
dynamics as well as of time-evolution. We recall that any radiative region of 
a rotating star is pervaded by baroclinic flows that appear in the first plac e 
as a differential rotation and a meridian circulation (e.g. Rieutordll2006al) . 
These flows control some turbulence there and affect the whole star, playing 
an important role in the transport of elements and momentum. 

We may now wonder what would be the ide al model o f an isol ated rapidly 
rotating star. This question was addressed in iRieutordl (|2006bl) and we re- 
produce here his discussion, which is still fully relevant: 

"Such a model should describe the mean state of a star at any time of its life 
and especially the new quantity specific to these stars: angular momentum. 
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"Unlike a non-rotating star, which is a one-dimensional object (in a large- 
scale description) which needs only scalar fields ([we] forget magnetic fields), 
a rotating star is, at least, a two-dimensional object with, at least, one vector 
field in addition to all scalar fields. Hence, complexity increases not only by 
the multi-dimensional nature of the model but also by the number of physical 
quantities to be determined. This implies that the ideal model deals consis- 
tently with angular momentum and especially the losses and gains through 
stellar winds and accretion. Such a model should also take into account the 
baroclinicity of radiative zones and there, the anisotropic turbulence which 
appears through shear instabilities; it should also include a mean-field the- 
ory of convection to forecast Reynolds stresses and heat flux. Of course, 
observers would like to know the emissivity of the atmosphere as a function 
of latitude (if they use interferometry) or as a function of wavelength if they 
do spectroscopy. But if they do asteroseismology they surely wish to know 
the eigenspectrum of these objects. 

"The foregoing points show that progress in the understanding of rotating 
stars needs also some advances in the following questions of stellar physics: 

• How angular momentum is distributed in a star and how is it input or 
output with what consequences ? 

• The immediately following question concerns the nature of the Reynolds 
stresses in the convective and radiative zones. 

• Then, what is the baroclinic state of the radiative regions ? 

• Similarly, the atmosphere is in a baroclinic state and cannot be at rest: 
how strong are the differential rotation and the meridional currents? Does 
the atmosphere develop strong azimuthal winds streaming around the 
star like Jupiter's winds? 

• Gravity darkening can be so efficient that equatorial regions are cool 
enough to develop convection; this raises the question of the latitude 
dependence of emis sivity of the atmosph ere beyond the von Zeipel model 
(see the attempt of Lovekin et al. 20061 ) . 



[we] did not mention magnetic fields. Clearly, they multiply the number 
of problems and first steps should ignore them if possible." 



3 Setting the problem 

Many of the foregoing questions can be answered by steady rotating stars, 
which do not evolve, neither dynamically by losing mass and angular momen- 
tum, nor by gravitationally contracting, nor by nuclear evolution. As a first 
step towards the general models, we concentrate on this simpler problem. 
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3.1 Equations for a steady rotating star 

We consider a lonely rotating star in a steady state. The star is governed by 
the following equations for macroscopic quantities: 



( A<(> = AirGp 
pT\r ■ Vs = -divF + e* 

p(2fl» A v + v Vv) = - VP - pV(<j> - \Qls 2 ) + F v 
div(pv) = 0. 



(1) 



The first equation is Poisson's equation which relates the mass distribution 
(p is the density) and the gravitational potential <j> (G is the gravitation con- 
stant). The second equation is the heat equation, which relates the advection 
of entropy s by the velocity field v when nuclear heat sources e* are present, 
with the heat flux F. The third equation is the momentum equation written 
in a frame rotating at angular velocity ri». P is the pressure field and the 
viscous force. Finally the last equation ensures mass conservation. 

These equations need to be completed by those describing the microphysics 
and the transport properties of stellar matter. We propose to describe the 
heat flux with 

F = - Xr VT-^^V S 
Km 

where \r is the radiative conductivity and Xturb a turbulent heat conductivity. 
In this expression the second term is assumed to represent the convective heat 
flux, which is supposed to be controlled by the entropy gradient. TZm = 1Z/A4 
where 1Z is the ideal gas constant and M. the mean molecular mass of the 
fluid. 

A realistic modelling of the viscous force would be derived from the turbu- 
lent Reynolds stresses, however before reaching this long term goal we assume 
the simple model of a constant dynamical viscosity p,, namely: 

F„ ^(ZW + ivdivv). 

The microphysics is completed by the three expressions 

(P=P(p,T) 
K = n(p,T) (2) 
£* = £*(/0,T) 

which respectively give the equation of state, the opacities and the nuclear 
heat generation. We recall that in radiative diffusive equilibrium, heat con- 
ductivity is related to opacity by 
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16crT 3 



Xr 



3np 

where a is the Stefan-Boltzmann constant. 



3.2 Boundary conditions — Angular momentum 
condition 

The previous system of partial differential equations needs to be completed 
by boundary conditions which affect the gravitational potential, the velocity 
field, the pressure and temperature. 

At the star centre these conditions are simply that all functions should be 
regular. At the surface things are more difficult. First because the surface of 
a star is ill-defined. To simplify, we shall assume that the surface is an isobar 
or an equipotential. If the properties of the model are independent of the 
chosen surface, then we may say that the choice of the bounding surface is 
not crucial. On this surface we have to: 

• match the gravitational potential to that in the vacuum, which vanishes 
at infinity; 

• impose stress-free conditions on the velocity field, namely that 

v • n = and ([cr]n) A n = 

where [a] is the stress tensor. The fluid is assumed to not flow outside the 
star, and to feel no horizontal stress (n is the outer normal of the star); 

• specify the pressure on the surface. A simple choice is P s = § = , where k 
and g are the averaged opacity and gravity on this surface; 

• impose the temperature boundary conditions. To simplify we assume that 
the star radiates locally as a black body. Therefore we ask 

n • VT + T/Lt = 

where Lt is specifying the scale of temperature variations near the sur- 
face. We note that if the diffusion approximation is used then Lt = 
I6/(3/5k), k being the opacity. 

We note that this problem, as formulated by Eqs. ([TJ, is not fully con- 
strained because the total angular momentum is not specified. We have spec- 
ified the rotation rate of the frame where the solution is computed but this is 
not specifying the rotation rate of matter, which is the combination of both 
v and the solid body rotation of the frame. To remove this degeneracy we 
may either specify the total angular momentum of the matter or specify the 
equatorial velocity of the star. In the first case we may write 
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r sin 9pu v dV = 

(V) 

and 

v v (r = R,6 = n/2) = 

in the second case. The first condition just specifies that all the angular 
momentum of the star is in the solid body rotation of the frame, while the 
second just says that the frame is rotating at the equatorial angular velocity. 



3.3 Scalings 

As formulated by ([T|), the problem is that of the steady flow of a self- 
gravitating compressible gas subject to some heat sources. As all fluid me- 
chanics problem we first need to choose the various relevant scales that feature 
the physical quantities. These scales are used to rewrite the set of equations 
with dimensionless variables. We converged to the following set: 



Length scale: equatorial radius R 

Time scale (R 2 /n AI T £ }^ 2 _ 

Velocity scale V = t/TZmT c 

Density scale: central density p c 

Temperature scale: central temperature T c 

Pressure scale 1ZmP c T c 

Entropy scale TZm 

Potential scale AnGR 2 p c 



This scaling uses the sound travel time as the time scale, and as a conse- 
quence the velocity field is scaled by a central sound velocity. Accordingly, 
the potential scale would be V 2 , but we prefer to use AttGR 2 p c , which comes 
from the free-fall time scale. 



3-4 Dimensionless equations and numbers 



Using the foregoing scaling leads to the following dimensionless equations 

1 



p (2He z A u + (u • V)u) = - Vp - pV(A p ip - -Q z s 2 ) 



E(Au+ -|Vdivu) 



pTu ■ Vs = div (xVT + XtTVs) + e 
div pu = 
{Aip = P 



(3) 
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with the numbers 

n = -2£=, e = a„ = ^ gr2 <- 



and the dimensionless functions 

Xr(p,T) _ Re*(p,T) 



a-kGr? Pc ' a n M PcRsFR^Tc p c (n M T c )W 

Although the critical angular velocity can only be computed after the com- 
pletion of the calculation, it is useful to define the pseudo-critical angular 
velocity f2 c — ^/AirGpc, 

^2 ^ 



and introduce other dimensionless functions, namely 

e(0) R 2 s»{p c ,T c ) 



s = e/s(0), X = X/X(0), A 7 



X(0) T cXr (p c ,T c )- 



3.5 Global parameters, p c ,T c ,R of the star model 

Let us suppose that we know the mass of the star M, how can the foregoing 
equations be used to determine the stellar quantities, especially p Cl T c and 
R? 

If we have solved the problem, we observe that A p is known, which gives 
a first relation between p c , T c and R 2 . In the same manner, the central value 
of e gives a relation between p c and T c : 

Re*{p c ,T c ) 



P c{n M T c f/ 2 - 

Thus, taking into account that the mass of the star expresses as a function 
of R and p c , namely 

M = p c R 3 ( pdV, 

J(V) 

the expressions of A p and £(0) completely determine the parameters of the 
stars, i.e. the radius, the central temperature and central density. 
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3.6 Numerical method 



The resolution of the set of equations ([3]) imposes some numerical challenges. 
First, the shape of a rotating star is spheroidal and not known a priori. We 
thus use a mapping of coordinates adapted to this geometry and modify 
it during the calculation. Secondly, as the problem is two-dimensional, its 
size quickly grows with the resolution, thus imposing the use of an efficient 
numerical method that should be, at the same time, precise enough to deal 
with velocity fields. Finally, as the problem is non-linear, we solve it using an 
iterative procedure. 



3.6.1 Computational domain and mapping of coordinates 



Following iBonazzola et al.l (|I998j ). we use a mapping of the coordinates r(£, 9) 
adapted to the geometry of the star. The new spheroidal coordinates (C,9), 
are defined in such a way that £ = 1 represents the surface of the star, 8 
being the colatitude. Thus doing, the surface of the star is a surface of coor- 
dinate, a property that very much simplifies the implementation of boundary 
conditions. 

The star can be subdivided in several domains, and in each domain the 
relation between spherical and spheroidal coordinates is given by the general 
expression 

r = a,( + Ai{C)[Ri+i{0) - airji+i] + Bi(()[Ri(6) - ana] m < C < m+u (4) 

where -Bi(C) = 1 — A'(C)- Here Ri(9) and Ri + i(9) represent the internal and 
external boundary of the domain respectively, and a,i, Ai(() are chosen to 
satisfy some properties at the boundaries between the different domains. In 
particular, we want 

r(C = m,0) = Ri(9) and r(( = m+1 , 9) = R l+1 {8), 

then it should be that Ai(rji) = and Ai = I. The value of the param- 
eters di should be such that r is monotonically increasing with ^. 

A nice property of this mapping is that it reduces to the spherical co- 
ordinates near the centre. The use of a spherical harmonic expansion of the 
unknowns is therefore possible, and the asymptotic properties of the solutions 
near the centre are well-known. The decomposition of the stellar domain into 
multi-subdomains improves the efficiency of the spectral methods to be used, 
especially in dealing with discontinuities (interface between a convective and 
radiative region) and in dealing with the large pressure and density variations. 

As we have already mentioned, at the stellar surface, the gravitational 
potential must match the vacuum solution vanishing at infinity. There is 
no easy way of imposing this condition on a surface with arbitrary shape. 
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Fig. 1 The computational domain: the star at the centre is divided into sub-domains 
(layers) that ease the computations. It is surrounded by a vacuum domain limited by a 
sphere which allows an easy implementation of the boundary conditions on the gravitational 
potential. 



To circumvent this difficulty, wc follow iBonazzola et al. I (Il998h who use an 



additional domain such that 1 < ( < 2 and whose outer boundary is a sphere 
as shown in Fig. Q] On this sphere bounding the computational domain, the 
gravitational potential can meet simple boundary conditions for each of its 
spherical harmonics, namely 

dfa | (1+1)0, _ Q 
dr r 

which ensure the matching with a field vanishing at infinity. Fig. [T] shows the 
full computational domain along with the key surfaces of the mapping. 

To write down the equations (j3|), the operators should be written with 
the spheroidal coordinates. Let first us recall that the relation between these 
coordinates and the classical spherical coordinates is 



r = r(C,0'), e = 9', <p = <p'. 



(5) 
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In order to manipulate these new coordinates, and especially the vectors, we 
use the natural basis defined by 



dr dr dr 



Making explicit these definitions, we can express the covariant and contravari- 
ant vectors, namely: 

E c = r c e r , E e = r g e r + re g , = r sin 9e v 



r^ rr^ r r sin 9 

as functions of the usual unit vectors associated with spherical coordinates 
(e r , eg, e v ); we set: 

dr dr 

Passing from the spherical to the spheroidal coordinates may be done using 
the following relations: 

d£ = ]_d£ dl^dl ndl d£ = df_ 

dr r c d( ' d9 d9' r c d( ' dip dip' ' U 

These expressions give the metric tensor whose covariant (gij = Ej • Ej) 
and contravariant (g y = E* • E J ) components are: 

9CC = r o 9C0 = r C r e, gee = r 2 + r 2 e , g w = r 2 sin 2 6> (7) 

9 a = T ^, 9 C6 = -^, 9 e6 = \, 9^ = ^~ (8) 
r z r£ r A r^ r z r 2 sin 6 

and g Cv = g 8ip = g Cip = g 6ip = 0. At r = 0, g et \g vv are singular. The metric 
determinant is such that 

yj~\g\ — r 2 r^ sin9 and \e t: > k \ = ^ 



The volume element is given by 

dV = y/\g\d(d6dip = r 2 r c sm.6dQd6dip = r 2 r c dC,dQ. 

As an example, the Poisson equation for the gravitational potential takes the 
form 
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Atp = g cc 



P, 



d 2 j> 



2r e d 2 xp 



1 



r 2 rQ 8(89 r 2 



rr C J 
1 8 



ree + r g cot ( 
r 2 r^ 



sin 



1 8 2 ip~ 
^8^_ 



dtp 

W 



where 



— and r i e ~ ~§c§9' 



no 



3.6.2 Numerical method 

As we have mentioned earlier, due to the 2D nature of the problem, the size of 
the calculation will increase very quickly with the number of grid points. This 
motivates the use of a high order method that can achieve high precision with 
low resolution grids, optimizing computation time and memory requirements. 

An additional difficulty comes from the fact that the overall structure of 
the star (profile of pressure, temperature, etc.) and its dynamics (rotation, 
meridional circulation) must be computed simultaneously. This requires a 
great precision and the ability to calculate higher order derivatives of some 
variables. 

Sp ectral methods are s pecially well-suited for this kind of problem ( GrandclementJ 



2006t Canuto et al. 20061 ). These methods expand the solutions on a basis of 
orthogonal functions. The approximation of the solution using n basis func- 
tions will be 

n— l 

<j>W(x) = Y,<*iPi(x)- (9) 

1=0 

The basis functions Pi(x) used are usually a set of orthogonal polynomials, 
such as the Legendre or Chebyshev polynomials. In our case it can be shown 
that if (p(x) is infinitely differentiable, the approximate function <f^ n > (x) will 
converge to the exact solution <p(x) faster than any power of the grid resolu- 
tion h. This expresses the fact that, in spectral me thods, the error decreases 
exponentially with the number of basis functions n ( Fornber3|l9 98) . The i-th 
derivative of a function is approximated by 



[axi =2>drf- (10) 

v ' 1=0 



Due to the non-linearity of the equations that we have to solve, we have 
to calculate the product of two variables of the model in an efficient way. 
Unfortunately, multiplication cannot be easily performed using the spectral 
coefficients, as it involves a convolution in the transformed space. This can 
be solved by using an alternate approach of spectral methods, with the same 
properties, called pseudospectral collocation methods. In a pseudospectral 



14 



Michel Rieutord and Francisco Espinosa Lara 



method the variables are not represented by their spectral coefficients, but 
by their values at certain points Xi called the collocation points. Then, all the 
calculations can be performed in the real space. 

The basis functions are orthogonal against some scalar product 
(Pi,P m ) — Si m . Then, we can get the spectral coefficients ai by 

ai = (<t>(x),Pt(x)). 

The scalar product (usually a weighted integral over the interval) can be 
calculated using the gaussian quadrature formula associated with the family 
of orthogonal polynomials Pi 

n-l 

ai =^2w j Pi(x j )(f)(x j ), 

1=0 

where Xj and Wj are the nodes and weights of the gaussian quadrature. Note 
that the collocation points. Then, the first derivative at the collocation 

points can be obtained as 

71— 1 /n-l \ 
n—l sn—1 \ 
n-l 

= Y,D ij <j>(x j ) 

3=0 

where Dij is the differentiation matrix. 

We use this procedure in the radial and latitudinal directions to transform 
the original system of non-linear partial differential equations in a system 
of non-linear algebraic equations. For the latitudinal direction, we use Leg- 
endre polynomials as the basis functions, while for the radial direction we 
use Chebyshev polynomials associated with the Gauss-Lobatto collocation 
points. This grid includes the points at the extrema in order to deal with 
boundary conditions. 

The nice convergence properties of spectral and pseudospcctral methods 
are only valid for smooth functions which are infinitely differcntiable. How- 
ever, it is known that inside a star there will be some discontinuities, as for 
example, at the boundary between a convective core and a radiative enve- 
lope. This difficulty is solved by using a multi-domain approach, in a way 
that the variables are continuous and differentiable within each domain but 
not necessarily at the boundaries between different domains. 
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3.6.3 Iterative procedure 

The system of algebraic equations resulting from the discretization of the 
problem is non-linear and is solved using an iterative method. For that we 
either use the well-known Newton's method or a relaxation method. For 
Newton's method we write the problem in the form 

F(x) = 0, (11) 

where the vector function F represents the equations that we want to solve 
and x is the vector containing all the independent variables of the problem 
(pressure, temperature, . . . ) including the shape of the surface which is not 
known a priori. The equations are linearized using the Jacobian matrix of 
F(x) defined as 

<5F(x) = J(x)<5x. (12) 

Then the correction to the solution in the fe-th iteration will be calculated 
solving the linear system 

J(x (fc) )<5x (fc) = -<5F(x (fe) ) (13) 

and we set x (fe+1) = x 1 ^ + 6x^ k K 

With an appropriate initial approximation y^°\ Newton's method has 
quadratic convergence. In practice, a rotating stellar model can be calculated 
in approximately 10 iterations starting with the corresponding non-rotating 
model. 



3.6.4 The case of flow fields 

The computation of the flow field is certainly the most delicate part of the 
solution. Its computation needs to circumvent two difficulties: first the flow 
faces very large variations of the density (typically eight orders of magnitude) 
and second the low viscosity of the stellar fluid. Indeed, this latter complica- 
tion implies that the flows need to be computed within the asymptotic regime 
of low Ekman numbers while the zero viscosity solution is degenerate for the 
linear part of the velocity op erator (any geostrophic flow may be added to a 
solution see Rieutordl 2006a! for detailed e xplanations). We shall present th is 



rather technical point in a separate work ( Espinosa Lara &: Rieutordl 2012 ) 



3. 7 Tests of the results 



We checked the results with two global tests: the virial theorem and the global 
balance of energy. 
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3.7.1 The virial test 

Let us first present the virial test. For this, we recall that the equations of a 
steady flow in a rotating frame are: 

2ft A pu + pu ■ Vu = -pVcf) + pfi 2 se s + Div[cr] 
div pu = 

with the boundary conditions on the velocity field u • n = and n A [er]n = 0. 
The virial equality is obtained by integrating the scalar product of r with the 
momentum equation over the fluid's volume. Hence, we have to evaluate the 
following integrals: 

• The z-component of the relative angular momentum 

/ r-2ftApudV = -2ft- / r A pudV = —2QL Z . 

J(V) J(V) 



• The gravitational energy 



r- P V4>dV= I I P 4>d 3 r = W. 



(V) 2 J{v) 



• The kinetic energy due to bulk rotation as measured by the frame rotation 
/ r • P Q 2 ae s dV = [ ps 2 dV = W 2 , 

J(V) J(V) 



l(V) An 

where / is the moment of inertia. 
The relative kinetic energy 



r • pu • Vu dV = - / pu 2 dV = -2T re i. 

(V) J(V) 



• The stress integral 

/ r • Div[er] dV = / ?*j<7jj dSj — / a„ dV. 

J(V) J(S) J(V) 

Stellar gas is assumed to be a newtonian fluid without volume viscosity, 
hence the stress tensor is 

(Tij — f-lCij pfiij 

where [c] is the shear tensor (cij = diVj + djVi — 2(dkVk)Sij /3), so that 
era = -3p- 

Thus the virial equality may be written 
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2T rc i + ml + W + 3P + I s + 2f2,L z = 
or, with non-dimensional quantities 

2T rol + IQ 2 + A P W + 3P + I s + 2QL Z = 
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(14) 



where I s is the surface integral that appears in the stress integral. In the 
case of a steady configuration like the one under consideration, the surface 
integral is estimated as follows: 



/ yLTidj dSj - \ 
J(S) J(; 



pr-dS 



(S) 



with 



n = E c / ||E C 



r • n 



r cVg' 



a' 



dS = r^r 2 + r 2 e dQ. 



Using nondimensional quantities, this leads to 



fin . ( 



2 r e (I du r dug u g . 

2d r u r - - divu- - [-— + — --) ) -p(9) 



r 3 dQ. 



3.7.2 The energy test 

Another test of internal coherence of the solutions is provided by the energy 
balance between sources and losses. The integral of the entropy equation over 
the fluid's volume gives: 

/ P Tu-VsdV= [ {xVT + xtTVs) ■ dS + [ edV. 

J(V) J(S) J(V) 

In the case of radiative envelopes at zero Prandtl number, this equation can 
be simplified to 



(xVT)-dS+ / sdV = 0. 

(S) J(V) 



4 Some results 



The foregoing algorithm has been used to compute models of rotating stars in 
a steady state. For these models we use an analytic expression for the energy 
generation rate, namely: 



£* (p, T, X, Z) = e Q {X, Z)p 2 T- 2 ^ cxp (,4/T 1 / 3 ) 
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Fig. 2 Differential rotation in a 5 Mq stellar model rotating at Q = 0.7i?x- 



as in lEspinosa Lara fc Rieutordl (|2007| ) . It is completed by the use of OPAL 



tables for the computation of the opacity and the derivation of the density 
from the equation of st ate (X = 0.7 and Z = 0.02 with solar composition of 
ICxrevesse fc Noelslll993h . 



We computed a few models of stars with a convective core, assumed to be 
an isentropic region, and with a radiative envelope. Presently, no convective 
envelope can be included in the models, which are therefore limited to stars 
with masses above 1.5 M Q . 

Numerical difficulties come both from the high density and pressure con- 
trast between centre and from high rotation rates. At the time of writing 
(April 2012) the most extreme models deal with 

Psmf/Pcontro = 10~ 14 and S7/S7 K = 0.9. 

The latter rotation corresponds to a flattening of 0.3. For these models, the 
spatial resolution uses a spherical harmonics series truncated at L = 64 and 
N = 400 collocation points on the radial grid, which are distributed over 8 
domains. 

As far as the velocity field is concerned, we recall that the baroclinic flows 
that pervade the radiative region are computed in the asymptotic limit of 
vanishing Ekman numbers. Thus, the Prandtl number is also set to zero 
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Fig. 3 The surface rotation as a function 
of colatitude. The stellar model is the same 
as that of Fig. [2] 



lell 




Fig. 4 Streamlines of the meridian circula- 
tion associated with the differential rotation 
shown in Fig. [2] Solid lines show counter- 
clockwise circulation and dashed lines clock- 
wise one. 



and heat advection by meridional flows is neglected. Fig. [5] illustrates the 
differential rotation that is forced b y the baroclinic torque. As also ob- 
served in previous (simpler ) mode ls of lEspinosa Lara fc Rieutordl (|2007l ) and 
Rieutord fc Espinosa Lara! ( 20091) . the core is rotating faster than the enve- 



lope. The differential rotation is cylindrical in the isentropic core (as required 
by the Taylor-Proudman theorem) , and almost shellular in the inner part of 
the radiative envelope. 

The meridional circulation shown in Figgis dominated by the streamlines 
of the Stewartson layers lying along the tangent cylinder of the convective 
core. This flow pattern should not be taken at face value since the interior 
flows are computed without viscosity. The balance of forces in the Stewartson 
layer cannot be ensured and therefore the flow pattern depends on the grid 
resolution. 

In Fig. [5j we provide a view of the squared Brunt- Valsala frequency. This 
shows the anisotropy of the buoyancy force which, especially in the outer 
layers, influences the g ravito-inertial modes that a re part of the oscillation 
spectrum of such stars ( Dintrans fe Rieutord 2000h . 

Besides the dynamics, we have also investigated the thermal structure of 
rotating stars and, mor e specifically, the distribut i on of the heat flux as a 
function of latitude. In lEspinosa Lara fc Rieutordl (|201ll ). 2D models have 
been used to validate a very simple model of the latitudinal variations of the 
flux, which depends on a single parameter Q/Qk- Such models are based 
on the idea that within an envelope the heat flux F satisfies div F = and 
is almost anti-parallel to the local effective gravity g e . Since the mass dis- 
tribution inside massive stars is concentrated the Roche model can be used, 
leading to latitudinal flux variations that only depend on the rotation rate. 
This simple model has been successfully compared with the very few obser- 
vations that are available and with complete two-dimensional models (see 
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Fig. 5 The square of the Brunt- Vaisala frequency distribution in a meridian plane. 



Espinosa Lara fc Rieutordl 120111) . As shown in Fig. [5] the two models nicely 
match and notably differ from the model of von Zeipel which predicts that 
T e ff oc gl^ 4 as a consequence of neglecting of the baroclinicity of the config- 
uration. 

In Tables Q] and [2] we show the physical parameters obtained from the 
calculations of a series of stellar models with masses between 3 and 20 Mq 
and rotation rates up to 90% of the breakup angular velocity. The models 
have been calculated using L max = 64 and 400 radial points distributed over 
8 domains. The virial and energy tests give an idea of the quality of the 
solutions. As we can see, the virial test values are very small, below a few 
10~ 10 , as a result of the high precision of spectral methods. The energy test is 
also good enough (under around 10 -5 ), the larger values being a consequence 
of the use of tabulated opacities which have limited precision. 



5 Conclusions 

The computation of two-dimensional stellar models is a real numerical chal- 
lenge which has been taken up by the ESTER project. The key difficulty is 
to derive simultaneously the bulk structure of the star and the mean velocity 
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Fig. 6 Latitudinal flux variation (gravity darkening) at the surface of a rapidly rotating 
stellar model of 5Mq . The solid line represents th e full two-dimensional model, the dashed 
line the simplest model described in text and in llEsp inosa Lara & Ricutord 201ll) . while 
the dotted line shows the prediction of von Zeipel hypothesis. 



fields that pervade it. These large-scale flows come from either the baroclin- 
icity of the radiative zones or from the anisotropic Reynolds stresses in the 
convective zones (although some Reynolds stresses might also be expected in 
the radiative zones if sh ear instabil ities can generate some small-scale tur- 
bulence as suggested bv lZahdll992h . At the moment, the ESTER code can 
produce dynamically self-consistent models, which include the background 
flows but no Reynolds stresses, for stars with masses larger than 1.5M© and 
rotation rates less than 90% the break-up. The possibilities of the code have 
not yet been fully explored however. 

As far as steady solutions are concerned, the main challenges are to enable 
the modelling of the lower-mass stars (solar type) with an outer convection 
zone, and to take into account the effects of viscosity and Reynolds stresses 
in the bulk of the stars. These latter effects might indeed be crucial to the 
transport of elements. 

The next important step will be to deal with time evolution; this should 
include: 



• dynamical evolution during PMS phase with a slow gravitational con- 
traction inducing spin-up of the star; 
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Table 1 Fundamental parameters for a series of rotating stellar models". 



M(M ) 17/17*. 


R(R ) 


e 


c 


n S ot(d) 




D eq (km 


./b) L(L q ) 


T cff (10 3 K) 


log 9e 


3.0 


0.0 


1.97 


0. 


00 






0.0 


81.2 


12.36 


4.33 


3.0 


0.3 


1.96(p) 
2.05(e) 


0. 


04 


0.70(p) 
0.65(e) 


0.54 


158.6 


80.0 


12.50(p) 
11.97(e) 


4.33(p) 
4.25(e) 


3.0 


0.5 


1.95(p) 
2.19(e) 


0. 


11 


0.46(p) 
0.43(e) 


0.36 


255.5 


78.4 


12.69(p) 
11.31(e) 


4.34(p) 
4.11(c) 


3.0 


0.7 


1.94(p) 
2.42(e) 


0. 


20 


0.38(p) 
0.36(e) 


0.30 


340.4 


77.0 


12.84(p) 
10.36(e) 


4.34(p) 
3.86(e) 


3.0 


0.9 


1.93(p) 
2.74(e) 


0. 


29 


0.35(p) 
0.34(e) 


0.28 


411.5 


76.4 


12.92(p) 
8.91(c) 


4.34(p) 
3.32(e) 


5.0 


0.0 


2.62 


0. 


00 






0.0 


542.8 


17.23 


4.30 


5.0 


0.3 


2.60(p) 
2.72(c) 


0. 


04 


0.83(p) 
0.77(c) 


0.66 


177.7 


533.2 


17.44(p) 
16.69(c) 


4.31(p) 
4.23(e) 


5.0 


0.5 


2.58(p) 
2.91(c) 


0. 


11 


0.55(p) 
0.51(c) 


0.44 


286.5 


520.9 


17.70(p) 
15.76(e) 


4.31(p) 
4.09(e) 


5.0 


0.7 


2.57(p) 
3.21(e) 





20 


0.45(p) 
0.42(e) 


0.36 


381.9 


510.4 


17.92(p) 
14.43(e) 


4.32(p) 
3.83(e) 


5.0 


0.9 


2.56(p) 
3.63(e) 


0. 


30 


0.42(p) 
0.40(e) 


0.33 


461.6 


505.3 


18.02(p) 
12.40(e) 


4.32(p) 
3.30(e) 


10.0 


0.0 


3.87 


0. 


00 






0.0 


5733.6 


25.55 


4.26 


10.0 


0.3 


3.84(p) 
4.02(e) 


0. 


04 


l-07(p) 
0.98(e) 


0.86 


206.8 


5619.0 


25.85(p) 
24.73(e) 


4.27(p) 
4.19(e) 


10.0 


0.5 


3.81(p) 
4.29(e) 


0. 


11 


0.70(p) 
0.65(e) 


0.57 


333.6 


5469.6 


26.24(p) 
23.35(e) 


4.27(p) 
4.05(e) 


10.0 


0.7 


3.78(p) 
4.73(e) 


0. 


20 


0.57(p) 
0.54(e) 


0.46 


444.8 


5341.5 


26.57(p) 
21.35(e) 


4.28(p) 
3.80(e) 


10.0 


0.9 


3.76(p) 
5.37(e) 


0. 


30 


0.53(p) 
0.51(c) 


0.43 


536.7 


5279.3 


26.73(p) 
18.28(e) 


4.28(p) 
3.26(e) 


20.0 


0.0 


5.70 


0. 


00 






0.0 


43791.2 


35.00 


4.23 


20.0 


0.3 


5.66(p) 
5.91(c) 


0. 


04 


1.36(p) 
1.24(c) 


1.13 


241.0 


42921.1 


35.41(p) 
33.89(e) 


4.23(p) 
4.16(e) 


20.0 


0.5 


5.61(p) 
6.32(e) 


0. 


11 


0.89(p) 
0.82(e) 


0.74 


388.5 


41775.4 


35.94(p) 
31.97(c) 


4.24(p) 
4.01(e) 


20.0 


0.7 


5.57(p) 
7.00(e) 


0. 


20 


0.72(p) 
0.68(e) 


0.61 


516.9 


40790.9 


36.39(p) 
29.16(c) 


4.24(p) 
3.76(e) 


20.0 


0.9 


5.54(p) 
8.02(e) 


0. 


31 


0.67(p) 
0.65(e) 


0.56 


621.0 


40326.9 


36.60(p) 
24.82(c) 


4.25(p) 
3.22(e) 



a From left to right: Mass, equatorial angular velocity, radius, flattening, central rotation 
period, surface rotation period, equatorial velocity, luminosity, effective temperature, 
logarithm of effective gravity (cgs). The values for the solar parameters used in the tabic 
are M = 1.9891 -10 33 g, Rq = 6.95508 -10 10 cm and L Q = 3.8396 -10 33 crg/s. 



c Flattening e = 1 - ^ . 
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Table 2 Fundamental parameters for a series of rotating stellar models (cont.) a . 



lVJ.llVJ.f7 


i O lOis 
) / i K 


p c (dyn/cm 2 


) Pc\B> 


Tc(K) 


Ps/Pc 




Ps/pc 




Virial 

test b 




Energy 
test c 




o.U 


U.U 


1.34 


10 17 


4U.O 


2.43 


10 7 


9.61 


10" 


15 


2.03 


10" 


ii 


5.08 


■10" 


ii 


5.84 


10" 


7 


3.0 


0.3 


1.34 


10 17 


41.0 


2.43 


10 7 


9.92 


10" 


15 


2.07 


10" 


11 


8.01 


■10" 


11 


3.42 


10" 


7 


3.0 


0.5 


1.35 


10 17 


41.3 


2.42 


10 7 


1.03 


10" 


14 


2.09 


10" 


11 


2.03 


■10" 


10 


6.14 


10" 


6 


3.0 


0.7 


1.35 


10 17 


41.5 


2.42 


10 7 


1.07 


10" 


14 


2.08 


10" 


11 


2.67 


■10" 


10 


7.50 


10" 


6 


3.0 


0.9 


1.36 


10 17 


41.6 


2.42 


10 7 


1.09 


10" 


14 


1.93 


10" 


11 


2.31 


■10" 


10 


1.09 


10" 


5 


5.0 


0.0 


8.00 


10 16 


21.3 


2.76 


10 7 


3.86 


10" 


14 


6.22 


10" 


11 


1.20 


■10" 


10 


3.75 


10" 


6 


5.0 


0.3 


8.04 


10 16 


21.4 


2.75 


10 7 


3.97 


10" 


14 


6.30 


10" 


11 


8.41 


■10" 


11 


8.84 


10" 


6 


5.0 


0.5 


8.09 


10 16 


21.6 


2.75 


10 7 


4.11 


10" 


14 


6.34 


10" 


11 


7.51 


■10" 


11 


7.95 


10" 


6 


5.0 


0.7 


8.13 


10 16 


21.7 


2.75 


10 7 


4.23 


10" 


14 


6.21 


10" 


11 


8.01 


■10" 


11 


7.33 


10" 


6 


5.0 


0.9 


8.16 


10 16 


21.8 


2.74 


10 7 


4.29 


10" 


14 


5.58 


10" 


11 


1.06 


■10" 


10 


1.34 


10" 


5 


1 n n 
1U.U 


U.U 


4.31 


10 16 


n a a 
y.44 


3.20 


10 7 


1.53 


10" 


13 


1.79 


10" 


10 


8.96 


■10" 


12 


1.12 


10" 


fi 


10.0 


0.3 


4.33 


10 16 


9.50 


3.20 


10 7 


1.56 


10" 


13 


1.80 


10" 


10 


4.91 


■10" 




5.53 


10" 


7 


10.0 


0.5 


4.35 


10 16 


9.58 


3.19 


10 7 


1.60 


10" 


13 


1.78 


10" 


10 


5.18 


■10" 


11 


1.25 


10" 


6 


10.0 


0.7 


4.38 


10 16 


9.65 


3.19 


10 7 


1.64 


10" 


13 


1.69 


10" 


10 


2.59 


■10" 


11 


5.62 


10" 


7 


10.0 


0.9 


4.39 


10 16 


9.69 


3.18 


10 7 


1.66 


10" 


13 


1.34 


10" 


10 


4.54 


■10" 


11 


4.48 


10" 


7 


20.0 


0.0 


2.79 


10 16 


4.88 


3.62 


10 7 


3.26 


10" 


13 


2.41 


10" 


10 


1.54 


■10" 


11 


2.38 


10" 


fi 


20.0 


0.3 


2.80 


10" 


4.91 


3.61 


10 7 


3.31 


10" 


13 


2.36 


10" 


10 


7.25 


■10" 


11 


1.49 


10" 


6 


20.0 


0.5 


2.81 


10 16 


4.94 


3.61 


10 7 


3.39 


10" 


13 


2.23 


10" 


10 


1.79 


■10" 


11 


2.11 


10" 


7 


20.0 


0.7 


2.82 


10 16 


4.98 


3.60 


10 7 


3.46 


10" 


13 


1.93 


10" 


10 


7.16 


■10" 


11 


4.38 


10" 


8 


20.0 


0.9 


2.83 


10 16 


5.00 


3.60 


10 7 


3.49 


10" 


13 


1.17 


io- 


10 


4.09 


■ io- 


11 


1.09 


10" 


6 



a From left to right: Mass, equatorial angular velocity, central pressure, central density, 
central temperature, ratio between surface and central pressure, ratio between surface 
and central density, relative virial test and relative energy test. The values for the solar 
parameters used in the table are M Q = 1.9891 ■ 10 33 g, R e = 6.95508 ■ 10 10 cm and 
L© = 3.8396 -10 33 crg/s. 

b Virial test: (2T rcl + IQ^ + W + 3P + I s + 2f!±L z )/W. 
- Energy test: (/ (S) (£VT) ■ dS + / (v) edv) / f (v) edV 



• dynamical evolution during main sequence with stellar wind, which forces 
a spin down and associated mixing; 

• nuclear evolution. 

These steps need the right algorithm for temporal evolution, which is not 
known presently. Indeed, such an algorithm should be at the same time fast, 
stable and precise. We have managed to use spectral methods, which ensure 
rapidity and precision but the stability remains a challenge. A better under- 
standing of the properties of the discretized operators is certainly a key to 
improve the efficiency of the algorithms. 



21 



Michel Rieutord and Francisco Espinosa Lara 



As shown by the foregoing examples, some realistic models can now be 
computed for intermediate mass and massive stars. These models are steady 
and therefore the chemical composition must be given. To circumvent this 
constraint, one-dimensional models can be used to compute time evolution 
and the ensuing chemical composition. Then, the bulk relation between pres- 
sure and chemical composition shown by the ID-model can be inserted into 
the 2D-model. In this way, steady models can include some stellar evolution. 

The steady models described in this work are most relevant in the study of 
the oscillation properties of rotating stars. The interpretation of the frequency 
spectrum of such stars is indeed a challenging problem and an intense use of 
2D-models will be necessary to find out how to invert data coming from, for 
instance, 5-Scuti stars. 

Another use of these models is obviously the interpretation of interfero- 
mctric data collected on some nearby fast rotating stars (a Aql, a Cep, a Leo, 
/? Cas, etc.). Fast rotating stars have a surface brightness that strongly de- 
pends on latitude (gravity darkening) . Accurate models are crucial to extract 
the physical parameters contained in the intcrfcrometric data from such stars. 

Finally, steady two-dimensional models may also serve as proxies for the 
initial conditions of a collapsing massive star, although the final word will 
come from time-evolved models including mass-loss. 
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